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INTRODUCTION 


Magnetic  circuits  are  currently  used  in  projectile  fuzing  circuits. 
Specifically,  magnetic  setback  generators  are  used  to  power  fuzing  cir- 
cuits . Computer  analysis , using  magnetic  field  vector  quantities , can  be 
a valuable  tool  in  the  optimization  of  such  systems . A theoretical  basis  of 
the  finite  difference  method  used  to  calculate  the  magnetic  vector  potential 
in  a nonlinear  magnetic  circuit  is  presented  herein . 

A typical  magnetic  setback  generator  is  shown  in  figure  1 . This 
configuration  will  be  used  to  derive  results  that  will  be  generally  appli- 
cable . 


The  determination  of  the  magnetic  field  distribution  of  a circuit 
containing  iron,  current  densities,  and  permanent  magnets  is  not  amen- 
able to  linear  analysis.  The  nonlinear,  partial  differential  equations  and 
boundary  conditions  for  the  cross-sectional  region  of  a machine  can  be 
defined,  but  analytical  solutions  to  the  defined  problem  are  only  possible 
for  the  most  simple  cases.  The  difficulty  is  due  to: 

Material  inhomogeneity  - There  are  large  differences  in  the 
reluctivities  of  the  magnetic  material  (iron)  and  nonmagnetic  materials 
(air , copper)  . 

Geometrically  complex  media  interfaces 

Material  nonlinearity  - The  reluctivity  is  a function  of  the 
flux  density  in  the  magnetic  material . 

The  approach  presented  here  is  to  transform  the  system  equations 
into  a finite  difference  formulation,  defined  at  a finite  number  of  discrete 
points  within  the  circuit.  The  resulting  set  of  simultaneous  equations 
can  then  be  solved  by  an  iterative  procedure,  while  taking  into  account 
the  material's  nonlinearity. 

The  following  assumptions  are  made: 

- The  field  distribution  is  two-dimensional . 

- The  field  is  static  (stationary)  . 

- The  machine  is  infinitely  long  (in  the  Z direction)  and  end 
effects  can  be  neglected . 
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Figure  1.  Magnetic  setback  generator 
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- The  magnetization  characteristic  of  the  materials  is 
represented  by  a single  valued  function,  i.e.  hysteresis  is  neglected. 

- The  material  is  isotropic  . 

- The  eddy  current  losses  may  be  neglected  (see  appendix  A)  . 

The  assumptions  allow  the  performance  of  a two-dimensional  analysis. 
Of  course,  a three-dimensional  model  would  be  more  representative  of 
the  physical  problem  (at  the  cost  of  increased  model  complexity  and  com- 
putation time)  . Nevertheless,  the  model  presented  should  yield  some 
useful  results . 

The  following  notations  will  be  used  throughout  this  report. 

A - Magnetic  vector  potential 
B - Magnetic  flux  density  (magnetic  induction) 

H - Magnetic  field  intensity 

I - Current 

J - Current  density 

L - Path  segment 

M - Magnetization  vector 

S - Surface  area 

ax  - Unit  vector  in  the  X direction 

h - Mesh  cell  dimension 

p - Permeability 

v - Reluctivity 

p - Charge  density 
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- The  vector  operator  (V  = 


ax 


^ a 

a + — 

X 3Y 


a 

az  az 


THE  NONLINEAR  VECTOR  POTENTIAL  FORMULATION 
The  Basic  Equation 

Electromagnetic  field  problems  can  be  solved  through  solution  of 
Maxwell's  equations  with  appropriate  boundary  conditions.  Unfortunately, 
in  many  cases  the  mathematics  become  difficult  when  using  Maxwell's 
equations  directly . The  magnetic  vector  potential  (A)  is  used  in  this 
formulation  to  simplify  the  mathematics . 

The  partial  differential  equation  describing  the  vector  potential 
for  steady  magnetic  fields  in  two  dimensions  using  rectangular  coordinates 
is: 


_a 

ax 


(v  9A  } + _9  (v  9A} 

1 ax  1 3Y  1 3YJ 


-J 


(1) 


where  v,  A,  and  J are  the  reluctivity  (see  appendix  B)  , the  Z component 
of  the  vector  potential,  and  the  Z component  of  the  current  density,  respec- 
tively . 


Equation  (1)  will  now  be  derived.  The  magnetic  field  intensity 
(H)  is  defined  to  be 

H = i |=vB  (2) 

where  p is  the  permeability 

v is  reluctivity 
B is  the  magnetic  flux  density 

The  point  form  of  the  Maxwell-Ampere  CURL  equation  may  be  stated: 

CURL  H = V x H = J (3) 


4 


The  magnetic  vector  potential  (A)  is  defined  by: 

CURL  A = V x A = B (4) 

and  from  Coulomb's  convention 

DIV  A = V • A = 0 (5) 

Since  we  are  only  considering  direct  currents , for  which  the 
charge  density  (p)  is  constant,  the  continuity  equation  for  the  current 
density  vector  (J)  is  expressed  as: 


DI  V J = v • J = - ^ 
dt 


= 0 


(6) 


That  is,  the  current  or  charge  per  second,  diverging  from  a small  volume, 
per  unit  volume  is  equal  to  the  time  rate  of  decrease  of  charge  per  unit 
volume. 

Substitution  of  (2)  and  (4)  into  (3)  results  in: 

CURL  II  = CURL  (vB)  = CURL  (v  CURL  A) 
therefore,  CURL  (v  CURL  A ) = J m 


It  is  now  assumed  that  the  system  under  study  possesses  axial 
symmetry , thereby  allowing  a description  of  the  magnetic  field  in  two 
dimensions,  i.e.  a cross  section  of  the  system  (appendix  C)  . It  is  also 
assumed  that  there  is  no  magnetic  flux  in  the  Z direction,  i.e. 


e-  _ 9B 

Br,  - 0,  „„  - 0.  In  addition,  it  is  assumed  that  current  will  only  flow 


JX  = 0’JY=0'  JZ  = J 


Z " 3Z 

perpendicularly  through  the  cross  section;  i.e, 

Since  the  vector  potential  is  in  the  same  direction  as  the  current  density, 

= = ® ant*  A=  A^.  Equation  7 can  be  written  in  its  explicit  form  by 

definition  of  CURL 


B = CURL  (A) 
f9A„ 


9A, 


,9Y 


9Z 


ax  + 


9A 
1 9Z~ 


9A. 


9Y 

(8) 


where  ax>  ay^  az  are  unit  vectors  in  the  X,  Y,  Z directions. 


Therefore , 


f9AZ 

9Ay) 

9Y 

9Z  j 

,9AX 

9Az\ 

9Z 

dX 

1 

9Ay 

9A  \ 
X 

9X 

" 9Y  j 

= B. 


= B. 


= B. 


Therefore  J = CURL  (v  B ) = CURL  (v  CURL  A ) 
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3Y 


3Z 


X 


1 3_ 

I 3Y 


Ir:  A 


0Axn 


3X  Y 9Y 


3A„  3 A. 


3Z  [ \ 9Y  9Z 
9 


9X 


9Ax  9Azn 


9Z 


9X 


+ 

9(v  Bx) 

3 (v  B 

9Z 

3X 

4. 

~9(v  By) 

8(vB? 

r 

9X 

3Y 

\ * 

ldAx 

3Az 

9Z 

V 

9Z 

3X 

_ 9 

3Ay 

3Ax| 

9X 

V 

3X 

3Y  j 

9 

- 

[9Az 

9Ay  , 

9Y 

V 

3Y 

3Z  i 

a,, 


i 


X 


>aY 


(9) 


yields 


Substituting  the  values  A - A - 0 and  A - A into  equation  (9) 
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dZ 

V 

~ 3X  j 
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It  is  assumed  that  J = and  = Jy  =0. 


Therefore , J = — 

z ax 


J = — 

x az 


z 

V|  ax 


3Ar 


a_ 

3Y 


,9A 

I — - 

■ 3Y 


/ j 


J = _ 

y az 


ax 


/ 3A, 

\ir 


= 0 


= 0 


Therefore,  f-  v — +—  v — 

ax  3X  3Y  3Y 


= - J 


(1) 


where 


J = J,  a, 

Z Z 


A = Az  az 


Equation  (1)  is  the  nonlinear  equivalent  of  Poisson's  equation. 

For  air  regions  which  do  not  have  current  density,  equation  (1) 
reduces  to  Laplace's  equation. 


32  A 32  A 


ax * a y 


r = o 


(ii) 


For  air  regions  which  include  the  current  density  fields , 
equation  (1)  reduces  to  the  linear  Poisson  equation: 

32A  32A  -J 

ax2  +3Yr"^“  (12) 

where  is  the  reluctivity  of  free  space. 


For  the  iron  regions  of  the  magnetic  circuit,  it  will  be  assumed 
that  the  current  density  is  zero,  i.e. , no  eddy  currents.  This  is  realistic 
if  we  consider  the  iron  to  be  made  of  thin  insulated  laminations  lying 
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normal  to  the  current  density  vector  . In  these  regions,  the  reluc- 
tivity is  a function  of  magnetic  induction  (B)  so  that  equation  (1)  be- 
comes: 


9 . 3A  , 9 , 9A  , . 

ax  (v  ax  } ay  (vay  5 0 


(13) 


Permanent  Magnet  Sources 


In  order  to  include  permanet  magnets  in  the  vector  potential  form- 

ulation,  the  concept  of  a magnetization  vector  (NO  is  used.  Briefly,  a mag- 
netized body  behaves  as  though  there  were  an  internal  current  with  density 

(CURL  Ml  analagous  to  J (see  appendix  D)  . Under  the  previously  stated 

assumptions,  M is  in  the  X-Y  plane  and  perpendicular  to  J^. 

For  a system  including  permanent  magnet  sources , the  magnetic 
flux  density  (B)  is  defined  to  be 


B = |a  (H  + M) 

Rearranging  this  equation  produces: 

H = (-)  B - M=  v B-  M 
\l 


(14) 


(15) 


where  His  the  magnetic  field  intensity 
B is  the  magnetic  flux  density 
Mis  the  magnetization. 

Restating  the  point  form  of  the  Maxwell  Ampere  CURL  equation 
(3)  shows: 

CURL  H = V x H = J 

Substituting  equation  (15)  into  equation  (3)  results  in: 

CURL  [vB-M]  = J (16) 
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Since  Vx  (A  + B)  = VxA+VxB,  equation  (16)  can  be  arranged: 
CURL  (vB)  = J + CURL  (M)  (17) 


By  substituting  equation  (4)  into  equation  (17)  : 

CURL  (v  CURL  A)  = J + (CURL  M) 

Previously  it  was  found  that,  under  our  assumptions 

CURL  (v CURL  A)  - ,v  || , + (v|f  , 

Therefore,  by  substitution: 


9 9A  9 3A  -*■ 

55  tv  ax  1 + 57  <v  57 > * ’[J  + (CURL «)] 


Next,  CURL  Mis  explicitly  expressed  as: 
CURL  M = 


(18) 


(19) 


f^z 

9My  \ 

1 j 

a + 

1 x 1 

,3Mx 

8Mz' 

9My 

aMx ' 

9Y 

9Z 

iaz 

3X  i 

aY  + 

i9X 

3Y  j 

(20) 


For  the  circuit  under  consideration 


9M  3M 

Mz  = 0a"d—  =0,^=0 


Also,  considering  only  the  X-Y  plane 


9M  9M 

= 0 — — 

9Z  ’ 9Z 


= 0 


therefore,  by  substitution 


CURL  M = 


9My  9Mx 


9X 


9Y 


(21) 
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the  total  expression  becomes: 


8 

8X 


. 8A  . 
(V0^ 


, 8A  . 


-J  + 


8Y 


8My 

ax 


(22) 


For  air  regions  not  having  current  density,  equation  (22)  reduces 
to  Laplace's  equation: 


= 0 


For  air  regions,  including  current  density,  equation  (22)  reduces 
to  the  linear  Poisson  equation 

82 A 82 A - J 

9X^  0Y2  * vQ 


For  the  iron  regions,  it  is  assumed  that  the  current  density  and  the  mag 
netization  vector  are  both  zero.  Therefore,  equation  (22)  reduces  to 


0 0A  0_ 

ax  lv  ax  J + 8Y 


, 3A  ^ 
v aY 


= o 


For  the  permanent  magnet  regions  assume  that  the  current  density  vector 
is  zero,  and  therefore,  equation  22  reduces  to 


DISCRETIZATION 


Grid  System 

In  order  to  use  numerical  solution  methods,  the  nonlinear  partial 
differential  equation  is  transformed  into  a difference  equation.  In  order 
to  do  this,  the  magnetic  circuit  cross  section  is  subdivided  by  horizontal 
and  vertical  gridlines.  The  vector  potential  is  then  evaluated  at  the  inter- 
sections of  the  gridlines,  yielding  an  array  of  discrete  point  values.  The 
rectangular  area,  formed  by  two  adjacent  horizontal  gridlines  and  two 
adjacent  vertical  gridlines,  forms  a mesh.  Nonuniformly  sized  meshes 
are  used  to  account  for  irregular  boundaries . That  is , a fine  mesh  is 
used  where  accuracy  is  needed,  such  as  boundaries.  Conversely,  a 
coarse  mesh  is  used  where  reluctivity  variations  are  small,  thereby  min- 
imizing the  number  of  gridpoints  and  resultant  computations . 

The  grid  system  is  constructed  so  that  each  row  is  cut  by  the  same 
columns  and  each  column  is  cut  by  the  same  rows . All  nonorthogonal 
boundaries  are  approximated  by  stair-type  contours.  A basic  mesh, 
formed  by  the  grid  system , is  shown  in  figure  2 . 

Reluctivity 


The  vector  potential  is  defined  at  the  grid  points  of  figure  2 , and 
the  reluctivities  are  defined  to  be  constant  over  the  area  of  each  of  the 
four  meshes  meeting  at  a typical  mesh  point  (0) . 


In  order  to  determine  the  reluctivity  of  the  rectangles,  0184,  0251, 

0362,  0473,  the  average  flux  density  (magnetic  induction)  at  the  center 
of  each  rectangle  is  calculated.  The  X 5 Y components  of  the  flux  density 

at  the  center  of  rectangle  0184  are  calculated  as 


(24) 

(25) 


(see  appendix  E)  where  h4  is  the  distance  between  meshpoints  0 and  4, 
and  hx  is  the  distance  between  meshpoints  0 and  1 . The  vector  potentials 
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Figure  2.  The  basic  mesh 
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at  points  0,  1 , 8,  4 are  denoted  as  A„  , Ax  , Ag  , A4  , respectively.  The 
absolute  value  of  the  flux  density  will  be 


(26) 


This  value  of  | B | is  used  to  calculate  the  reluctivity  by  using  an  approx- 
imation to  the  nonlinear  magnetization  characteristic  (see  appendix  B)  . 


The  reluctivities  in  the  other  three  rectangles  having  gridpoint  0 
as  a joint  corner  are  determined  similarly . 


FINITE  DIFFERENCE  FORM  OF  THE  VECTOR  POTENTIAL 

In  this  section,  the  finite  difference  form  of  the  vector  potential 
will  be  derived  using  two  approaches  , the  point  form  and  the  integral 
form  of  the  Maxwell-Ampere  circuital  law . Although  different  forms  are 
obtained  in  each  case,  the  results  will  be  shown  to  be  equivalent. 

Point  Form 


In  order  to  obtain  the  finite  difference  form  of  equation  (1) 


d_ 

ax 


x 3 , 3A  , _ _ 

+ 3Y  V dY  J 


In  terms  of  the  basic  mesh,  the  differentials  will  be  replaced  with  differ- 
ences (refer  to  figure  3)  . 


Let 


(9A) 

axil. 


(27) 


dA  _ A°  - A3 
3XJ_3  h3 


(28) 
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9A  A*  A° 

( — ) 

9Y  o h2 


) 29) 


(2A)  = 

9Y  l 4 


A0  A4 

hT" 


(30) 


Also  let 

Vab 

be  the  weighted  average  of  Vj  and  along  a-b 

Vbc 

be  the  weighted  average  of  v and  v along  b-c 

Vcd 

be  the  weighted  average  of  vIn  and  along  c- 

Vda 

be  the  weighted  average  of  and  Vj  along  d-a 

Therefore , 


ax  (vax  J 


ab 


3A 


ax 


- v . 
1-0  cd 


3A 

9X|  o_3 


hj  + hj 


f Aj  ~ A0  ^ 

> 

1 

> 

U 

1 

_VaJ 

i h‘  i 

- vcd  h‘  J 

hx  + h3 


(31) 


d_  dA  . 
3Y  (V  3Y  } 


1 

3A| 

(aYjj.o  Vda 

3Y 

0-4 


h2  + h4 


/ A2  Aq  1 

I A0  A4  | 

. be  l h2  j 

Vdal  h4  /_ 

h2  + h4 


(32) 
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Figure  3.  Mesh  cell 
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By  substitution  equation  (1)  becomes: 


(2)  Vab  (Al  " Ao)  " (2)  vcd  (A°  " + (2)  vbc  (A*  " Ao) 

hi  (hx  + h3 ) h3  (hx  + h3 ) h2  (h2  + h4 ) 


_ (2)  vda  (A0  - A4) 
h4  (h2  + h4 ) 


(33) 


Let 


® 2 


«3 


a4 


2 

hj  (hj  + h3 ) 
2 

h2  (h2  + h4 ) 
2 

h3  (hi  + h3 ) 
2 

h4  (h2  + h4) 


(34) 

(35) 

(36) 

(37) 


Substituting  into  equation  (33) : 

vab  (Aj  - A0 ) + vbc  a2  (A2  - A0 ) + Vcd  a3  (A3  - A0) 

+ Vda  “4  (A<  “ Ao)  = "JT 


(38) 


Solving  for  A0 : 

Ao 


JT  * A-a-vab  * A’°‘vbc  * A»a»vCd  ^ A<a,Vda 

a*  vab  + °2  vbc  + °3  vcd  + °4  vda 

(39) 
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The  average  reluctivities  can  be  written  as: 


Let 


V„v.  = vT  (t-  - ■ . ) + vTT  ( 


ab  I h2  + h4  II  h2  + h4 


) 


(40) 


= vTT  (r — — r— ) + vTTT  ( 


be  II  vhx  + h3  ' III  vhx  + h3 


(41) 


h4 


V = V ( ) + V ( 

cd  III  1 h2  + h4  } IV  1 h2  + h4 


(42) 


Vda  VIV  ( hx  + h3  } + VI  ( hx  + h3 


(43) 


where 


Vj  is  the  reluctivity  in  region  I (a  constant) 
is  the  reluctivity  in  region  II  (a  constant) 
is  the  reluctivity  in  region  III  (a  constant) 
is  the  reluctivity  in  region  IV  (a  constant) 
hi 


Pi  = 


hx  + h3 


(44) 


P2  = 


h2  + h4 


(45) 


Ps  = 


hx  + h3 


(46) 


P4  = 


h2  + h4 


(47) 
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The  expression  for  A0  becomes 


A0 


JT  + Al0l  (VjP4  + vnp2)  + A2a2  (vnPl  + vmp3) 

+ A3a3  (vmp2  + v^p,)  + A4  a4  ( y wp3  + vTpt) 

“1  (viP<  + V„P2 ) + a2  (VjjPi  + vmps ) + a3  (vmp2  + v^p,) 
+ a<(vlVp3  + v^) 


(48) 


If  the  constant  current  density  in  the  mesh  cells  is  J . J . J J 
then  J^,  may  be  written:  * ® ^ 


^1  b4  hj  h2 

T I ( h2  + h3 ) (h2  + h4 ) II  (h3  + h3 ) (h2  + h4 ) 


h2  h3 


+ J 


h3  h4 


+ J. 


Ill  (hx  + h3)  (h2  h4)  iv  (h2  + h3)  (h2  + h4) 


= J,  PlP<  + J[J  PlPi  + Jjj,  PlPl  + Jjv  P,P, 


(49) 


Next,  permanent  magnets  are  included  in  the  formulation.  To  do  this,  the 
concept  of  the  magnetization  vector  (M)  , which  is  explained  in  appendix  D 
is  used.  Figure  4 represents  the  basic  mesh  and  includes  the  magnetization 

vector  (M)  for  each  cell.  The  magnetization  vector  is  assumed  to  be  a 
constant  value  throughout  the  mesh  cell  and , for  the  given  assumptions , 
it  has  only  X and  Y components . 


In  the  previously  derived  expression  (equation  (22) ) which  included 
permanent  magnets , : 


9_ 

9X 


, 9A 

V 9X  9Y  (V  9Y  J 


-J 


+ 


9Y 


9my 

9X 
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Figure  4.  Mesh  cell  with  magnetization  vector  (M) 
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The  only  term  not  evaluated  was: 


9MX  9My 


9Y  9X 


Let 


Y 


ths  weighted  average  of  the  Y component  along  ab 
Y 

(M)cd  be  the  weighted  average  of  the  Y component  along  cd 

x 

(M)bc  be  the  weighted  average  of  the  X component  along  be 

X 

^^da  be  t^le  weighted  average  of  the  X component  along  da 
Therefore  at  point  0 


3Mx|  <M>bc  - (M> 


X 

da 


9X  ^ q h4  + h2 


(50) 


2 


(51) 


2 


By  definition 


(52) 


(53) 
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(54) 


(M) 


Y Y 

Y Mj  h4  + M h2 
ab  h2  + h4 


(M) 


Mra  h » + h< 


cd  h2  + h4 

Therefore  by  substitution 


~Mn  hl 

, x _ - 
+ Mm  h» 

M; % h3  + M*h/ 

aMx 

hi 

+ h3 

hi  + h3 

9Y 


h4  + h. 


(55) 


= 2I 

[Ki  - mXi)  pi 

+ l 

Mm  - O P>] 

h2  + h4 

Y Y 1 

Mj  h4  + Mn  h2 

Mm  h2  + Miv  h* 

mY 

h2  + h4 

h2  + h4 

ax 

1 h!  + h3  | 

1 2 / 

. 2[Ky-  Q< 

(hx  + h3 ) 


(56) 


(57) 
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Therefore,  the  complete  expression  is: 

Aiai  <viP4  + vnp2)  + A2a2  (Vjjp!  + vmp3)  + A3a3  (vinp2  + vJvp4) 
+ A4a4  (v^Ps  + Vjpj ) + JT  + Mt 
° ~ °i  (vjp4  + vnp2 ) + a2  (vnp!  + vmp3 ) + a3  (vmp2  + Viyp4 ) 

+ (v^Pa  + v pj  ) 

(58) 

where , 

2 

Ql  hj  (hx  + h8) 

2 

“2  h2  (h2  + h4) 

2 

as  h3(hj  + h3) 

2 

04  h4  (h2  + h4 ) 
hl 

P1  ~ h[  + hs 
h2 

p2  = h2  + h4 
h 3 

Ps  ~ T 7~ 

hi  + h3 

h4 

^ h2  + h4 
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[(MjY-  MY  )p4  + (MY  - MYn)p2] 


- [(M*  - M^Pj  + (M*n  - M*  )p3] 


(h2  + h4 ) 


Integral  Form 

Ampere's  circuital  law  may  be  used  to  derive  the  finite  difference 
form  of  the  vector  potential . Ampere's  circuital  law  may  be  stated:  The 
line  integral  of  the  magnetic  field  intensity  (H)  about  any  closed  path  is 
exactly  equal  to  the  current  enclosed  by  that  path  i.e. 


(59) 


Positive  current  is  defined  as  following  in  the  direction  of  advance  of  a 
right-handed  screw  turned  in  the  direction  in  which  the  closed  path  is 
traversed. 


Since 


B 


CURL  A 


(60) 


and  H e vB 


(61) 


(62) 


Therefore , H = v B = v CURL  A 
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Substituting  into  equation  (59) : 


-/ 


I = d)  v CURL  A • d L 


(63) 


Next,  apply  Ampere's  law  along  the  path  a-b-c-d-a  as  shown  in  figure  3, 
Therefore , 


/ 


b c 

v CURL  A • dL  = / v CURL  A • dL  + / v,  CURL  A • dL 
a ab  b bc 


+ / vcd  CURL  A • dL  + / v CURL  A • dL 
c d a 


(64) 


The  weighted  average  reluctance  along  each  of  the  path  segments  is: 


vIh*  + vuh2 
ab  h2  + h4 


(65) 


vnh‘  * vmh« 

be  hx  + h3 


(66) 


; _ vmh2  + Vivh4 

cd  h2  + h4 


(67) 


Viyh,  + V h, 

) = 

da  hx  + h3 


(68) 


Since  the  dot  product  between  two  orthogonal  unit  vectors  is  zero , the 
following  products  are  zero: 
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Therefore,  after  replacing  the  differentials  with  differences,  the  path 

segments  become: 

h2  + h4 

A L = 

ab  2 

h l + h3 


AL,  = 
be 


A L 


2 

h2  + hit 


cd  2 

hi  + h3 


A L 


da 


and 


CURL  A = 


[aA z' 

* J-8Az| 

3Y 

aY  1 ax  J 

(See  appendix  E . ) 


Therefore , 


/ v , CURL  A -dL  = 
a ab  ab 

a a 


/>  r-^z  - i 

f Vab  [ ax  aY  j * dLab 


= V 


ab 


r - aaz  i 

AX  I ALab 


vjh4  + vnh2  * 

> 

> 

o 

i 

r 

- h2  + h4 

1 

•C 

h2  + h4 


(Aj  Ae ) 

hi 


V,h4  vnh, 


(69) 


25 


c c 

/ v CURL  A • dL  = / V. 

- DC  DC  , bC 

b b 


[ 8Az  - 1 

-k- 

[AAzl 

. 3Y  ax. 

^^c  Vbc 

. ay 

AL. 


be 


vIIhl  + vIIIhs 

A2  “ A0 

/hi  + h,l 

- 

hi  + h3 

— _ 

h2 

.1  2 i 

(A2  - A0) 


vnhi  vmh3 

— i — + — 


(70) 


/ v , CURL  A • dL  = / v 

cd  cd  cd 

c c 


a±z_  - 
ax  aY 


dL  =v  . 
cd  cd 


AAzl 

sr  Kn 


vIIIh*  + vIVh< 

i 

r-  — 

Aq  “ A3 

- (h2  + h4) 

h2  + h4 

L 1 -j 

»3  - 

2 

(A0  - A,) 


v — h2 


v„,h4 


III"4  JJV 


(71) 
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/ v , CURL  A 
, da 
a 


dL,  = / v, 

da  , da 
d 


9A„  1 f AA_  ] 

Z -*■  ! Z 

a I • dL  = v 

9Y  X1  ’ 


viyh3  + Vjh!  1 
hj  + h3  j 


A0  - A4 


da  da  \ AY  j 
hx  + h3 


AL 


da  i 

j 


Aq  A4 

hT“ 


Vrvh3  vThi 


h4 


2 2 , 


(72) 


Assuming  that  the  current  density  (J)  in  each  region  of  the  mesh 
is  constant,  the  total  current  enclosed  by  rectangle  a-b-c-d-a  in  figure  3 is 


I0  = //  J • dS  = J_  (AS)  = 


Jjhxh,  Jnhxh2  Jmh2h3  JIV*13  h4 


4 (Jihlh<  + Jnhlh2  + Jmh2h3  + J!vh3h^ 


(73) 


Combining  previous  results 

I0  =J>  H • 


dL  = • 


1 

(Ax  - A0) 

r 

vjh4  Vjjhj 

(A*  - A0 ) 

VIIhl  , VIIIh3 

hi 

2 ' 2 

h2 

2 ' 2 

L_  . j 

r~ 

r- 

(A0  - A3) 


vinhg  + viyh< 


(A0  - A4 ) v h3  v hj 

— — +-K 


L 


J 


(74) 


Solving  for  the  vector  potential  at  point  0: 

4 


An  — 


I0+  2 A a! 
i=l  1 


E a! 
i=l  1 


(75) 
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where 


1°  ' 4 (Jihlh4  + Jnhlh2  + Jmh2h3  + 


a = 


Vjh4  + vnh2 

2hq 


ai  = 


vnhi  + vmh3 

2 h. 


as 


vmh2  + viyh< 

2h, 


ai  = 


vIVh3  + vlh* 


2 hA 


To  include  permanent  magnets  in  the  formulation,  the  concept  of 
the  magnetization  vector  is  used. 

The  magnetic  induction  (B)  in  a system  including  air , iron , currents , 
and  permanent  magnets  is  described  by 

B = p (H+  M)  (76) 


which,  rearranged  is 

1 -!»  — ^ 1 

H = — (B)  - M =vB  - M 

H- 

Substituting  equation  (77)  into  equation  (59) 
/(vB  -M  ) • dL  = I 


(77) 


(78) 
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or 


^ v B • dL  - ^ 


M • dL  = I 


(79) 


Rearranged  r ^ ^ r ^ 

vB*dL  = I+®M  • dL 


Since 


therefore 


B = CURL  A 


• / 


(v  CURL  A)  • dL  = I + 4>  M • dL 


/ 


(80) 

(81) 

(82) 


/ 

The  only  term  not  previously  evaluated  is^/i  M • dL  . Using 
figure  4 as  a basis , 


0 


M*dL  = / M • dL  + / M • dL,  + / M • dL 

_ ab  be  cd 

a b c 


+ / M • dL 


da 


(83) 


Note,  that  M in  the  two-dimensional  case  is  located  in  the  X-Y  plane,  i.e. , 

v- 

Since  the  dot  product  between  two  orthogonal  unit  vectors  is  zero; 
the  following  products  are  zero: 

M , • L = M l • L,  = M • L , = M , • L , =0 

ab  ab  be  dc  cd  cd  da  da 


(84) 
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The  weighted  average  of  M along  each  path  segment  is: 


Y Y 

„ Mi  h,+  M„  h„ 

m'Y  1 4 11  2 

Mab  = 


Y 

Mbc  = 


Y 

McS 


h + h 

2 

4 

mii  »i  * 

MIII  h3 

h + 

h 

1 

3 

my  h 

+ M Y h 

III  2 

IV  4 

h 

+ h 

2 

4 

X , 

X , 

M,„  h.  + 

Nf  h, 

IV  3 

I 1 

A V 

M X = h,+  h, 

da  13 


(85) 


(86) 


(87) 


(88) 


Therefore, 


b 

f 

a 


• dtab  = 

Y 

M h + 
I •» 

Y . 

M b 

II  2 

h„  + h 

2 4 

+ 

eg 

rC 

h4 

CM 

— — 

(89) 


/ 

b 


d 

f 

c 


MX 

II 

h 

1 

x . 

+ M h 

III  3 

h 

+ h„ 

- 

l 

3 

V 

h 

+ MY  h 

III 

2 

IV 

h 

+ h 

2 

4 

h3> 


2 


(h  + h ) 
2 4 


(90) 


(91) 
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/ M,X  • dL 
d da 


da 


„IVX  h3  * MX  h,  ' 

hi  + h3 

hl  * h3 

CM 

(92) 


and. 


M • dL  = 1/2 


M?  - mii K 


Y Y l 

M - M | h 
II  III  2 


MX-M  X 
IV  III 


\ * i>4  - mwK 


(93) 


Substituting  this  result  into  equation  (75),  the  final  expression 
for  the  vector  potential  becomes 


I Ai  + I0  + M0 
A - i=l 

Ao 


4 

E a: 

i=l  1 


(94) 


where  Vi\  + vIIll2  vII^l  + VIII^3 


•I 

2hl 

' “2  = 

2h. 

2 

. VIII 

+ 

CM 

X 

vIV^4  » 

_ a 

vIVh3  + vIhl 

“3 

CO 

X 

CM 

a4  “ 

2h4 

io 

= 1/4 

(>J  jh^h 

4 + 

Jnhi] 

i2  + 

JIIIh2h3  + JIV 

Mn= 

r 

= 1/2  S 

(mx  - 

Mx  ! 

lh  + 

my 

% 

0 

l 

\ 1 

u\ 

1 1 

II 

III)  2 

+ |m 

lx  • 

- Mx 

h ■ 

► |Mv  - mV  lh 

IV 

III 

3 

[ I ivj  4 
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The  formulation  derived  using  the  point  form  of  the  .Maxwell  Ampere 
law  is  equivalent  to  the  equation  just  derived. 

From  previously  defined  variables, 


1 2 1 

i!v  / **  u v i 

f h 2 \ 

[ h2+h4 ) 

|W9j 

= 2(vl\  +VIIh2) 
h1(h1+h3)  (h2+h4) 


(h1+h3)(h2+h4) 


vjh4  + Vjjh2 j 

2h~i  I 


(h1+h3)  (h2+h4) 


Cap 


(95) 


Similarly, 


“2  Cvnsi  + viiiB3)  * 


(h1+h3)  (h2+h4) 


Ca2) 


(96) 


a.  (vtttB„  + vt„B,.)  “ 


(a:) 


3 ''VIIIP2  "iv’V  " (h1+h3)(h2+h4)  v~3 


(97) 


a4  (vIVB3  + Vjf^)  = (h^+h^)  (h2+h4) 


(98) 


32 


Substituting  into  equation  (58) , 


An  = 


i &ia'i  + i0  + 
i=l 


(hj+hg)  (h2+b4)  ) 


M7 


4 

z of 
i=i  1 


(99) 


Next,  the  magnetization  term  is  evaluated. 


(hi+h3)  (h2+hit)  j (hi+h3)  (h2+h4) 


•j  Mr  - 


2[{Mj  - Mjv]  84  + (Mil  - Mux)  ^2] 

(VhJ 


-2  [|M*t  - ♦ jM^n  - M*v}63]  ■ 


<VV 


= (l/2)[(h1+h3)(h2+h4)] 


Y Y 
M - M 
I IV 


! Y Y 
4 l II  "HI 


6 + |M„  - M 


6 (h+h  ) 


2j  2 4 


K - Miii ♦ 


X X 

m - m ig 


III  IV / 3 j * 1 3 


h t-h 


(h  +h  ) (h  +h  ) 
v 1 3JV  2 

(1/2) 

f 

if 

Y Y I h4  j mY  mY  1 

h2  ] 

liM 

i “ Miv]  Ch2+K)  + |Mn  " miii 

<VV  J 

(h+h  ) 

2 4 


.X  ..X 


+ [{MI  ' MIl)(h,*h3)  * |MlV  MIII|  (h1+h3J 


<Vh3> 
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= 1/2 


MX  - MX  ]h  + 
I III  1 


MY  -my 
II 


+ (Mi  - MivJh- 


MX  - MX 


IV 


h 

III / 3 


= M. 


Cioo) 


Therefore,  equations  C58)  and  C94)  are  equivalent  i.e. 
4 

E + lo  + 


E a' 
i=l  1 


^lal  Cvj84+Vjj82)  + A2a2  Cvj j8 j+Vj j j63)  + A3a3  Cvjjj62+Vjy64) 

+ (Viv^3+Vl^l)  + *^T  + 

°1  ^VI^4+VH^2^  + a2  ^VII^1+VIH^3^  + a3  ^VIII^2+VIV^‘t^ 

+ CI4  CVjy83+Vj8j) 
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BOUNDARY  CONDITIONS 


In  order  to  solve  the  vector  potential  equations,  it  is  necessary 
to  determine  boundary  conditions. 

Using  figure  5 boundary  conditions  that  are  consistant  with  the 
assumptions  can  be  determined.  The  various  boundary  types  can  be 
described  as: 

Type  A - Line  A-B  is  an  axis  of  symmetry,  and,  in  this  case, 
it  is  also  a center  line  of  the  magnetic  circuit.  Therefore,  A-B 
represents  a line  of  flux,  i.e.  an  equipotential  where 

A-B  = constant  potential  (A) 

Type  B - The  equi-potential  A-B  has  its  return  path  B-C, 

C-D,  D-E,  E-F,  F-A 

Type  C - It  is  assumed  that  the  magnetic  induction  outside 
the  contours  of  the  machine  is  negligible. 

Type  D - This  type  of  boundary  exists  at  air- iron  interfaces. 
At  these^boundaries,  the  tangential  component  of  the  magnet  intensity 
vector  (H)  and  the  normal  component  of  the  magnetic  vector  (IT)  must 
be  continuous.  For  a boundary  that  is  parallel  to  the  X direction, 
these  conditions  are  satisfied  when: 


and 


,3A, 

Vl3YJ  iron 


,3  A. 

voCgy) 


air 


(3A)  _ r_3A-j 

3X  iron  “ 'ax  air 


At  air-iron  boundaries  parallel  to  the  Y-axis,  these  conditions  will 
be  satisfied  when: 


and 


/ 3A->  _ f 3 A-. 

v 3X^  iron  “ v°  ^alr  air 


(— ) . = (9A) 

3Y  iron  3Y  air 
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Boundaries  a -Z  and  b-c  are  of  this  type. 


Type  E - At  air-iron  interface  corners,  such  as  point  Z 
and  point  c in  figure  5,  all  four  conditions  of  Type  D must  be  satis- 
fied simultaneously. 

Type  F - This  type  boundary  is  at  the  air-current  field 
region  interface,  such  as  £-c.  It  is  assumed  that  the  current 
carrying  material  has  the  same  permeability  as  air.  At  the.  boundary, 
the  tangential  component  of  H and  the  normal  component  of  B must  be 
continuous.  Therefore,  at  the  boundary 


(— ) 
hXJ 

current 

C— ) 

^dXJ  air 

region 

(—■) 

'•3Y-' 

current 

region 

'•3Y-'  air 

Type  G - This  type  of  boundary  is  at  the  iron-current  field 
region,  such  as  c-d,  d-e,  h-k,  or  k-£.  Since  the  reluctance  of  the 
current  carrying  region  is  the  same  as  in  air,  the  boundaries  are 
the  same  as  Type  D boundaries. 

Type  H - This  boundary  is  at  the  permanent  magnet- current 
carrying  region  interface,  represented  by  line  e-h.  Since,  in 
figure  5,  this  boundary  is  parallel  to  the  Y axis,  the  following 
conditions  must  be  true  at  the  boundary: 


(.BA.  - (M.) 

^dYJ  magnet  V3Y' 


current 

region 


fd A, 

' dXJ  magnet 


^0 


(— ) 
{‘dXJ 


current 

region 
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Figure  5.  Magnetic  circuit  with  boundaries  noted 
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Type  K - This  boundary  is  at  the  permanent  magnet- iron 
interface,  i.e.  lines  e-f  and  g-h.  Since  the  boundaries  are  parallel 
to  the  X-axis,  the  following  conditions  must  hold: 


v ( _ c 

3Y  iron  V dY^  magnet 

(M.)  . = rM.) 

9X  iron  V9X^  magnet 


SOLUTION  PROCEDURES 

The  vector  potential  (Aq)  at  a grid  point  not  on  an  outer 
boundary  is  given  by  equation  (94),  evaluated  after  the  appropriate 
reluctivities  and  a and  (3  constants  have  been  substituted.  Thus,  N 
simultaneous  equations  are  obtained  for  the  potentials  at  N internal 
grid  points.  The  solution  of  the  set  of  N simultaneous  equations 
yields  the  numerical  value  of  the  vector  potential  at  the  N grid 
points. 

The  solution  of  the  N simultaneous  equations  can  be  obtained  by 
an  iterative  procedure  via  computer.  This  is  necessary,  because  the 
coefficients  in  equation  (94)  are  functions  of  the  reluctivities  that 
are  in  turn  functions  of  the  CURL  of  the  vector  potential  calculated 
by  equation  (94)  . Using  the  relaxation  method  requires  a two-step 
iterative  process . In  the  first  step,  the  coefficients  in  equation 
(94)  are  assumed  to  be  constant  and  the  vector  potential  is  calculated. 
In  the  second  step,  the  coefficients,  which  are  functions  of  the  re- 
luctivity, are  recalculated  from  the  numerical  values  of  the  vector 
potential  of  step  one.  The  new  values  are  then  used  to  calculate  the 
vector  potential.  This  process  is  repeated  until  the  vector  potentials 
have  converged  to  values  that  are  considered  close  enough  to  final 
values . 
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In  general,  the  solution  procedures  may  be  summarized 

1 . Input  of  the  permanent  magnet  data,  the  gridsystem, 
and  the  nonlinear  relationships  v (reluctivity )=  f(B)  and 

M (magnetization)  = f(B). 

2.  Calculation  of  the  coefficients  of  the  difference 
equations  (initially,  coefficients  are  assigned  constant  values) . 

3.  Solution  of  the  difference  equations  with  point  intera- 
tion,  line  iteration,  elimination  (i£  v and  M are^  constant)  or  with 
the  Newton-Raphson  method  (if  v and  M depend  on  B) . 

4.  Acceleration  of  convergence  if  the  difference  equations 
are  solved  by  iteration. 

5.  Check  of  convergence:  Steps  2-5  are  repeated  until  a 
specified  convergence  criterion  is  met. 

6.  Calculation  of  characteristic  values  of  the  magnetic 
field,  such  as  the  magnetic  induction,  the  vector  equi-potential 
lines  (lines  of  flux),  etc. 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  finite  difference  approach  to  analysis  of  nonlinear  electro 
magnetic  circuits  is  a viable  technique.  One  possible  formulation 
and  its  theoretical  basis  are  presented  in  this  report.  In  addition, 
system  equations  corresponding  to  those  utilized  in  the  existing 
CLYDE  computer  program  (appendix  F)  are  shown. 

It  is  recommended  that  the  investigation  be  continued.  Speci- 
fically, it  should  be  determined  if  the  CLYDE  program  can  be  modified 
to  solve  nonlinear  electromagnetic  circuit  problems,  or  if  a new 
computer  program  must  be  developed. 
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APPENDIX  A 


EDDY  CURRENT  LOSSES 


In  the  two-dimensional  model,  infinite  length  in  the  Z direction  is 
assumed  as  well  as  homogeneity  in  the  Z direction.  If  no  additional 
assumptions  are  made  a loss  of  energy  in  the  iron  regions  due  to  eddy 
currents  (circulating  in  the  X-Z  plane)  must  be  taken  into  account.  The 
introduction  of  laminations  reduces  this  energy  loss,  but  also  reduces 
the  effective  cross-sectional  area  of  the  magnetic  material. 


EDDY 


Therefore,  the  insulation  thickness  between  laminations  is  assumed 
to  be  zero.  Then,  the  only  effect  will  be  to  eliminate  eddy  currents. 
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APPENDIX  B 


APPROXIMATION  FOR  RELUCTIVITY  (v) 


In  order  to  utilize  the  finite  difference  method  for  calculating  the 
vector  potential,  the  nonlinear  magnetization  curves  of  the  materials 
under  consideration  must  be  represented. 

The  magnetization  curve  is  the  well-known  hysteresis  loop  shown 
in  figure  B-1 . This  curve  is  nonlinear  and  multivalued.  Since,  for  most 
materials,  this  hysteresis  loop  is  rather  thin,  it  is  assumed  that  the  mag- 
netization curve  is  monotonic  and  single  valued.  By  using  average 
values  of  H for  each  value  of  B (or  average  values  of  B for  each  value  of 
H),  average  curves  that  closely  follow  the  so-called  virgin  magnetization 
characteristic  can  be  obtained.  In  addition,  it  is  assumed  that  the  mag- 
netization curves  pass  through  the  origin. 


Therefore, 


HAVE 


H + H0 
u £ 


or  B 


AVE 


B 

u 


+ B 


i 
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A plot  of  the  average  magnetization  curve  is  shown  in  figure  B-2 . 


In  order  to  minimize  the  calculation  time  for  the  vector  potential, 
a piecewise  linear  approximation  is  used  for  the  average 
magnetization  curve.  As  shown  in  a figure  B-3,  it  is  assumed  that 
beyond  a saturation  magnetic  induction  (Bs) , the  permeability  is  that  of 
free  space.  The  distance  between  Bs  and  the  origin  is  subdivided  into 
n equal  parts,  denoted  as  AB.  Next,  all  Hj's  are  read  and  stored.  The 
slope  between  all  Pj's  are  read  and  stored.  The  slope  between  all  Pj's 

and  P.  , 's  is  calculated  as: 

J+1 


m. 

J 


H.  , 

= _i±l 


H. 


AB 
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In  order  to  calculate  the  reluctivity  for  a flux  density  less  than 
in  the  interval  j to  j+1,  the  following  expression  is  used. 


v 


H 

B 


(B-B.)m.  + H. 
J J l 

B 


To  calculate  reluctivity  at  a point  B > B 


S: 


v 


Hs* 


<B-BS)V0 

B 


where  v^  is  the  reluctivity  of  free  space. 


s and 
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Figure  B-1.  Magnetization  curve  (hysteresis). 


Figure  B-2.  Average  magnetization  curve. 
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B (FLUX  DENSITY) 


Figure  B-3.  Piecewise  linear  approximation  of  the  magnetization  curve. 
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APPENDIX  C 


THE  LAPLACIAN  OPERATOR 


In  Cartesian  coordinates 


V2A  = + ifA_  + Y,  z) 

9X2  9Y2  3Z2 

In  two  dimensions  this  reduces  to 


V2A  = 


* /(x,  Y) 


In  cylindrical  coordinates 


v2a  = _1  i_(R  Ja|+  J_ 
R 9R\  9R I R2 


92A 

9Z2 


By  expanding  the  first  term: 


1 

R 


9A 
9R  . 


92A  J_ 
9R2  + R 


Therefore,  V2A  = + 1.  + / 9 2 A \ _ r,R 

9Z2  9R2  R 9R  R2!  9(J)2|  1 ’ 

In  two  dimensions  this  is: 


V2A  = d^A  + d2A  + 1 9A  _ f 21 

R,  Z 9Z2  9R2  R » U ’ J 


9A 

9R 

Z,  <f.) 
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The  relationship  between  coordinate  systems  is  defined  by  the 
following  representation: 


r=Vx2  + y2 


Z = Z 

4>  = Tan-1  "x" 

X=  R Cos  <f> 
Y = R Sin  <t> 

Z = Z 
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Therefore,  in  general. 


V2A  • 

Note: 
If  Z 

If  <f>  = 

If  <f>  = 


32A  + _9^A  + ^A  = _9^A  + i^A  + 1_  8*A  + ^A 

9X2  9Y2  9Z2  9R2  R 9R  R2  9<()2  9Z2 


= 0, 


v2a  = 


92A 

9X2 


+ 


92A 

9Y2 


+ _1  _3A  + £ ^A 
9R2  R 9R  R2  9(f)2 

0,  then  Y = 0 


V2A  = 


92A 

9Z2 


92A  1 9A  92A 

— 1 + — - + 

9R2  R 9R  9Z2 

90°,  then  X = 0 


V2A  = 


92A  + _9^A 

9Y2  + 9Z2 


92A  + 1 jA  + 9^A 

9R2  R 3R  9Z2 
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APPENDIX  D 


THE  MAGNETIZATION  VECTOR 

The  basic  source  of  a magnetic  field  is  electrical  current.  The 
magnetic  field,  due  to  electrical  current  flowing  in  a loop  and 
observed  at  large  distances  compared  with  the  dimensions  of  the  loop, 
is  called  a dipole  field.  The  small  current  loop  is  referred  to  as 
a magnetic  dipole.  The  magnetic  dipole  moment  (m)  is  defined  as: 


m = I S n 


where: 

I is  the  loop  current 

S is  the  area  of  the  loop 

. 

n is  a unit  vector  in  the  direction 
determined  by  the  right  hand  rule. 

Electrons  surrounding  the  nuclei  of  atoms  in  matter  have  spin 
and  orbital  motion.  In  effect,  the  electrons  form  current  loops. 

When  viewed  macroscopically,  these  current  loops  appear  to  be  a 
continuum  of  magnetic  dipoles. 

When  the  individual  dipole  moments  in  a material  are  randomly 
oriented,  their  net  effect  is  zero.  However,  if  there  is  a coherence 
to  the  orientation  there  will  be  a net  effect,  and  the  material  is 
said  to  be  magnetized. 

The  magnetization  vector  (M)  can  be  conceived  as  the  net  magnetic 
dipole  moment  per  unit  volume,  i.e.. 


i . E In 

M = lim  

AV+O  AV 

The  unit  for  magnetization  is  amperes  per  meter  and  is  the  same 
dimensionally  as  the  magnetic  field  intensity  (H) . The  vector  $ 
will  be  zero  if  the  dipole  moments  are  random  and  nonzero  if  there 
is  some  coherence  of  the  dipole  moments. 


» 
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A magnetized  body  Jjehaves  as  though  there  were  an  internal  current 
density  equal  to  CURL  M and  analagous  to  a current  density  vector  (J) . 
CURL  M is  usually  referred  to  as  the  atomic  or  molecular  current  density. 
That  is  because  the  field  contributions  of  CURL  M result  from  spin  or 
orbital  motions  that  are  localized  on  the  atomic  scale. 

The  general  definition  for  B in  terms  of  its  current  sources  is: 


IT  = y(TT  + J?)  . 
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APPENDIX  E 


X AND  Y COMPONENTS  OF  THE  MAGNETIC  INDUCTION 


B = CURL  A = VXA 


9AZ  3 A, 

3Y 


aY  + 

3Z  I A \ 3Z 


3AV  3A, 
3X 


V 


^ a„ 


3X  3Y 


Since  it  was  assumed  that  Ay  = A^  = 0 


3A_ 

\ 

-3A_ 

B = CURL  A = 

Z | 

Z 

3Y  j 

+ 

1 aX  + 

i 3X  , 

*Y 


Therefore, 


3A  3A 

BX  = and  By  = - — — 

3Y  3X 


As  an  approximation  the  average  magnetic  induction  in  the  center 
of  the  mesh  cell  can  he  calculated: 


3AZ  AAz  [a0  * A, 


A4  + Agl  Ai  - As  + Aq  - A4 


3X 

AY 

h4 

-3AZ 

-AAZ 

__ 

J Al  + Ae  1 J Ao  + A4  1 

l 2 / \ 2 n 

3X 

AX 

hl  J 

2h. 


Aj  - Ag  + Ag  - A4^ 


2h, 


Therefore 


* I' B I - yj Bx2  + By2 


-a8+Ao-a,A  2 
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Ai-Ao+Ag-A4 

2h, 
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APPENDIX  F 


CLYDE  SYSTEM  EQUATIONS 


Consideration  is  being  given  to  adapting  the  computer  program 
CLYDE1  to  handle  nonlinear  electromagnetic  circuits. 


The  formulation  of  the  Harmonic  operator  for  an  irregular  grid  in 
the  CLYDE  program  is: 


2A 

t + 

m ^ 

2B 

£ + 

2k 

h + 

2B 

b1(b1+b3) 

a 

b2Cb24b4* 

b 3(b  i+b  3) 

JO  T 

b 4 (b  2+b  4) 

2A  + 

b lb  3 


2B 

b 2b  4 


h = h2D 


A corresponding  equation  can  be  obtained  by  rearranging 


4 

A = E A a"  + I + M 
0 i-1  i i 0 0 

V 

Z a" 
i=l  i 


and  then 
h2V2A0 
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ltiplying  by  h2: 
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